#ifndef _MT3D_H
#define _MT3D_H

/*The followings bypass a bug that VsCode's IntelliSense reports incorrect errors when using eigen3 library*/
#if __INTELLISENSE__
#undef __ARM_NEON
#undef __ARM_NEON__
#endif

#include "gctl/core.h"
#include "gctl/io.h"

#include "Eigen/Dense"
#include "Eigen/Sparse"

#include "lcg/solver_eigen.h"
#include "lcg/preconditioner_eigen.h"
#include "omp.h"

#include "../utilities/glni.h"
#include "../utilities/efem_sf.h"
#include "../utilities/layerEarth_MT.h"


#define pi 3.1415926535897932384626433832795
#define _mu0 (4e-7*pi)
#define _eps0 8.8541878176e-12

typedef std::complex<double> complex_d;
typedef Eigen::SparseMatrix<complex_d, Eigen::RowMajor> spmat_cd; // 注意Eigen默认的稀疏矩阵排序为列优先
typedef Eigen::Triplet<complex_d> triplt_cd;
typedef Eigen::VectorXcd vector_cd;

static const complex_d _one = complex_d(1.0, 0.0);
static const complex_d _zero = complex_d(0.0, 0.0);
static const complex_d _1i  = complex_d(0.0, 1.0);

static int le_index[12] = {0, 1, 0, 2, 0, 3, 1, 2, 1, 3, 2, 3}; // local edge index

struct conductive_region
{
    int tag;
    double cndt, mu, epsl; // cndt (conductivity) epsl (epsilon)
    double ref_cndt; // background conductivity
};

struct mt_station
{
    gctl::vertex3dc *vert;
    std::vector<gctl::tetrahedron*> host;

    complex_d impedance[2][2]; // impedance tensor
    complex_d VMTF[2]; // vertical magnetic transfer function
    double phi[2][2]; // phase tensor

    //complex_d Zxy, Zyx;
};

struct tet_edge
{
    unsigned int id;
    gctl::vertex3dc *vert[2]; // Two ends of the edge
    std::vector<gctl::tetrahedron*> host; // Tetrahedrons that contain the edge

    tet_edge(){}
    tet_edge(gctl::vertex3dc *v1, gctl::vertex3dc *v2){set(v1, v2);}

    void set(gctl::vertex3dc *v1, gctl::vertex3dc *v2)
    {
        vert[0] = v1; vert[1] = v2;
        return;
    }

    bool operator ==(const tet_edge &e) // Overload the equality operator
    {
        if (e.vert[0] == vert[0] && e.vert[1] == vert[1]) return true;
        return false;
    }
};

struct edge_setup
{
    unsigned int uv_id[2];
    efem_sf::ordertype_e uv_type[2];
};

class mt3d : public CLCG_EIGEN_Solver // Import the solver class
{
public:
    mt3d(){}
    virtual ~mt3d(){}

    void run(int argc, char *argv[]); // Run the process

    void read_config(std::string filename); // Read configuration
    void read_model(std::string filename); // Read tetgen model and initialize arraries
    void save_model(std::string filename, bool debug = true); // Save Gmsh files
    void save_observation(std::string filename); // Save MT observations

    void get_element_conductivity(int tag_num, double &conduct, bool ref_cndt = false); // Get element conductivity by tag number
    void assemble_linear_system(); // Calculate A and B for the linear system Ax = B
    void apply_boundary_condition(); // Modify the linear system to enforce the boundary condition
    void solve(); // Solve the linear system
    void post_process(); // Calculate the MT observations

    void check_validity(); // Help function: Test to see if the A is symmetric

    // 偏微分方程的左端项
    static double _Aij_part1(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att, void *instance) // 辅助函数 类的成员函数的外部接口
    {
        return reinterpret_cast<mt3d*>(instance)->Aij_part1(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,att);
    }

    static double _Aij_part2(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att, void *instance)
    {
        return reinterpret_cast<mt3d*>(instance)->Aij_part2(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,att);
    }

    double Aij_part1(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att) // 偏微分方程的左端项1
    {
        edge_setup *et = static_cast<edge_setup*>(att);
        double ux, uy, uz, vx, vy, vz, dummy;

        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[0],efem_sf::Curl,et->uv_type[0],ux,uy,uz,dummy,dummy,dummy);
        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[1],efem_sf::Curl,et->uv_type[1],vx,vy,vz,dummy,dummy,dummy);
        return (ux*vx + uy*vy + uz*vz);
    }

    double Aij_part2(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att) // 偏微分方程的左端项2
    {
        edge_setup *et = static_cast<edge_setup*>(att);
        double ux, uy, uz, vx, vy, vz, dummy;

        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[0],efem_sf::Value,et->uv_type[0],ux,uy,uz,dummy,dummy,dummy);
        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[1],efem_sf::Value,et->uv_type[1],vx,vy,vz,dummy,dummy,dummy);
        return (ux*vx + uy*vy + uz*vz);
    }

    static std::complex<double> _Bi_TE(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att, void *instance) // 辅助函数 类的成员函数的外部接口
    {
        return reinterpret_cast<mt3d*>(instance)->Bi_TE(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,att);
    }

    std::complex<double> Bi_TE(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att)
    {
        edge_setup *et = static_cast<edge_setup*>(att);
        double ux, uy, uz, dummy;

        le_te_.calculate_EH(-1.0*z, 0.0);
        //le_.calculate_EH(-1.0*z, 0.0);

        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[0],efem_sf::Value,et->uv_type[0],ux,uy,uz,dummy,dummy,dummy);

        // le_所使用的坐标系为纵轴向下的右手直角坐标系
        // 此程序使用的是纵轴向上的右手直角坐标系
        return (ux*le_te_.get_Ey() + uy*le_te_.get_Ex() - uz*le_te_.get_Ez());
        //return (ux*le_.get_Ey() + uy*le_.get_Ex() - uz*le_.get_Ez());
    }

    static std::complex<double> _Bi_TM(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att, void *instance) // 辅助函数 类的成员函数的外部接口
    {
        return reinterpret_cast<mt3d*>(instance)->Bi_TM(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,att);
    }

    std::complex<double> Bi_TM(double x, double y, double z,
        double x1, double x2, double x3, double x4, 
        double y1, double y2, double y3, double y4, 
        double z1, double z2, double z3, double z4, void *att)
    {
        edge_setup *et = static_cast<edge_setup*>(att);
        double ux, uy, uz, dummy;

        le_tm_.calculate_EH(-1.0*z, 0.0);
        //le_.calculate_EH(-1.0*z, 0.0);

        elsf_.tetrahedron(x,y,z,x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4,
            et->uv_id[0],efem_sf::Value,et->uv_type[0],ux,uy,uz,dummy,dummy,dummy);

        // le_所使用的坐标系为纵轴向下的右手直角坐标系
        // 此程序使用的是纵轴向上的右手直角坐标系
        return (ux*le_tm_.get_Ey() + uy*le_tm_.get_Ex() - uz*le_tm_.get_Ez());
        //return (ux*le_.get_Ex() - uy*le_.get_Ey() - uz*le_.get_Ez());
    }

    // Realize the virtual function of the solver class
    virtual void AxProduct(const Eigen::VectorXcd &x, Eigen::VectorXcd &prod_Ax, 
        lcg_matrix_e layout, clcg_complex_e conjugate)
    {
        if (conjugate == Conjugate) prod_Ax = A_PDE_.conjugate() * x;
        else prod_Ax = A_PDE_ * x;
        return;
    }

    virtual void MxProduct(const Eigen::VectorXcd &x, Eigen::VectorXcd &prod_Mx, 
        lcg_matrix_e layout, clcg_complex_e conjugate)
    {
        prod_Mx = P_.cwiseProduct(x);
        //clcg_solve_lower_triangle(L_, x, P_);
		//clcg_solve_upper_triangle(LT_, P_, prod_Mx);
        return;
    }

protected:
    size_t node_num_, elem_num_, edge_num_, face_num_;
    gctl::array<gctl::vertex3dc> mesh_nodes_; // 网格中的顶点
    gctl::array<gctl::tetrahedron> mesh_elems_; // 网格中的四面体
    gctl::array<gctl::triangle> mesh_faces_; // 网格中的三角形面元（用于识别边界上的面与边）
    std::vector<tet_edge> mesh_edges_; // 网格中的边
    gctl::matrix<size_t> elem_edge_index_; // 四面体到边的列表
    gctl::matrix<bool> elem_edge_reversed_; // 四面体到边的索引记号，为真则表示对应边的索引在计算时需要反转
    gctl::array<int> node_markers_, elem_markers_, face_markers_, edge_makers_; // 元素的标签值
    gctl::array<bool> edge_on_boundary_;
    std::map<std::string, unsigned int> edge_index_map_;
    std::map<std::string, unsigned int>::iterator edge_index_iter_;

    glni::tetrahedron g3d_tet_;
    efem_sf::linear elsf_;
    spmat_cd A_PDE_, L_, LT_;
    vector_cd X_TE_, X_TM_, B_TE_, B_TM_, P_;

    std::string tetgen_file_; // 网格文件名
    double omega_; // 电磁波角频率
    double theta_in_; // 主场入射角
    int station_marker_; // 观测点标签值
    int bound_marker_; // 边界顶点标签值
    int air_marker_; // 空气层标签值
    std::vector<conductive_layer> primary_model_; // 主场模型 (导电率 相对磁导率 相对介电常数 底深度)
    gctl::array<conductive_region> secondary_model_; // 磁场模型 (标签值 导电率 主场导电率 相对磁导率 相对介电常数)

    size_t obs_num_;
    std::vector<mt_station> obs_site_; // 观测站

    layerEarthMT le_te_, le_tm_; // 层状介质的电磁场计算对象
    //layerEarthMT le_; // 层状介质的电磁场计算对象

    bool modeling_only_;
};

#endif // _MT3D_H

